% First set the initial parameters. 
pi=2*acos(0);                  % defines the constant pi using inverse cosine 
p=0;                           % start of computational domain 
q=1;                           % end of computational domain 
t=0;                           % start time 
Kx=1.0;                        % diffusion coefficient in x direction 
N=51;                          % number of grid points 
dx=(q-p)/(N-1);                % spatial step size 
x = p: dx : q;                 % vector of grid points 
u0=zeros(1,N);                 % vector of initial u values filled with zeros 
u0=sin(pi*x);                  % initial condition 
initialu=u0;                   % store initial u profile for comparison 
F=0.9;                         % safety factor 
dt=dx*dx/(2*Kx);               % von Neumann time step for stability 
dt=F*dt;                       % adjusted time step for safety 
ntimesteps=300;                % number of time steps 
Rx=Kx*(dt/(dx*dx));            % constant in scheme
% 
u=zeros(1,N+2); % define correct sized numerical solution array to include ghost values
% 
% Begin the time marching algorithm
disp('start time marching ...')
for counter=1:ntimesteps
 t=t+dt; % current time for outputted solution 
 u0=[0 u0 0]; % insert ghost values
 for i=2:N+1
 u(i)=u0(i) + Rx*(u0(i+1)-2*u0(i)+u0(i-1)); 
 end 
 u0=u(2:N+1); % copy solution to initial conditions for next iteration 
end % of time loop
disp('time marching ends, see graphics window for results')
exact=exp(-Kx*pi*pi*t)*sin(pi*x); % exact solution 
plot(x,initialu,x,exact,x,u(2:N+1),'k+') % plot of numerical and exact soln and initial profile
xlabel('x')
ylabel('concentration u') 
title('pure diffusion: initial profile (--) and numerical (o) and exact solutions')
